There are 13 questions to this activity. Save your answers in word document that you will hand in on Moodle using a .pdf extension. Keep your script file in your GitHub folder and make sure that all changes are pushed to GitHub. Code for all questions should be clearly commented with the question number. You will include a link to your script file as a part of your final question in this activity. The assignment is worth 50 points.
Glacial retreat will drastically affect the regional hydrology, vegetation, and the energy balance. For example, a surface with snow or a glacier reflects more incoming solar radiation off of the earth’s surface compared to bare rock and vegetation nearby. The loss of glaciers means that the land surface will experience even greater warming as more solar energy is absorbed. As glaciers retreat and overall snow and ice cover decreases in mountainous regions, vegetation can grow in newly formed bare surfaces. You can see in the picture above, there is grass growing in areas that were once covered in snow and ice. This growth of vegetation can drive additional change on the land surface and alter water, energy, and carbon cycling. A study by Carlson et. al. in 2017 demonstrated that less snow cover during the year drives increased vegetation in the French Alps. Another recent study by Karen Anderson et. al. in 2019 found that the Himalayas are experiencing increased vegetation growth as glaciers and permanent snow cover retreats. Read through these papers posted in the data folder to get an idea of the type of analysis and considerations for this type of research.
In this exercise, you’ll use data from the National Park Service that shows the footprint of 39 glaciers in 1966, 1998, 2005, and 2015 to show the change in glacier extent over recent decades. We’ll also look at a satellite derived measure of vegetation called the Normalized Difference Vegetation Index (NDVI). NDVI Values ranging from 0 to 1 indicate the relative amount of photosynthetic material on the ground. Values close to one indicate a higher level of vegetation and zero indicates no vegetation. The overall value will depend on characteristics of vegetation, and areas completely covered in vegetation do not necessarily correspond to an NDVI value of one. NDVI varies seasonally and typically peaks in the summer for many cold winter ecosystems. The maximum NDVI value indicates the highest amount of vegetation that year. Tracking changes in NDVI throughout time can indicate if there is increasing or decreasing amounts of vegetation.
You will use a series of spatial packages. Install the following packages:
#install.packages(c("raster","sp","rgdal","rgeos","plyr"))
You’ll notice there is a package plyr on the list. It is an older version of the package dplyr. We can’t use dplyr here because it overwrites too many spatial functions. If you were to load dplyr, you would get the following messages:
That is a lot of useful spatial functions that are overwritten by dplyr functions. You definitely want to avoid this when working with spatial data. Installation may take a few minutes. Next load all the packages to your environment.
library(raster)
library(sp)
library(rgdal)
library(rgeos)
library(plyr)
You will start by reading in the vector data of the glacier outlines. These glaciers have been digitized by the National Park Service. In shapefiles, all metadata is contained in the .xml file. You can open it in a web browser to find more information. The GNP_glaciers file contains all four shapefiles for each year. Read them in using the readOGR function.
#read in shapefiles
#readOGR in rgdal does this
g1966 <- readOGR("data\\GNPglaciers\\GNPglaciers_1966.shp", stringsAsFactors = T)
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\mloranty\Documents\GitHub\GEOG331_answers\activity6\data\GNPglaciers\GNPglaciers_1966.shp", layer: "GNPglaciers_1966"
## with 39 features
## It has 13 fields
## Integer64 fields read as strings: OBJECTID
g1998 <- readOGR("data\\GNPglaciers\\GNPglaciers_1998.shp", stringsAsFactors = T)
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\mloranty\Documents\GitHub\GEOG331_answers\activity6\data\GNPglaciers\GNPglaciers_1998.shp", layer: "GNPglaciers_1998"
## with 39 features
## It has 13 fields
## Integer64 fields read as strings: OBJECTID
g2005 <- readOGR("data\\GNPglaciers\\GNPglaciers_2005.shp", stringsAsFactors = T)
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\mloranty\Documents\GitHub\GEOG331_answers\activity6\data\GNPglaciers\GNPglaciers_2005.shp", layer: "GNPglaciers_2005"
## with 39 features
## It has 13 fields
## Integer64 fields read as strings: OBJECTID
g2015 <- readOGR("data\\GNPglaciers\\GNPglaciers_2015.shp", stringsAsFactors = T)
## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\mloranty\Documents\GitHub\GEOG331_answers\activity6\data\GNPglaciers\GNPglaciers_2015.shp", layer: "GNPglaciers_2015"
## with 39 features
## It has 21 fields
## Integer64 fields read as strings: OBJECTID Recno
Note there are 39 features (glaciers) for each shapefile and a table with information about each glacier that has 13 columns of data. Let’s investigate what the format of this data.
str(g2015)
#data stores all accompanying info/measurements for each spatial object
head(g2015@data)
## OBJECTID Recno Year Src1_0822 Src1_nadir Src3_nadir Src2_nadir Area2015
## 0 1 1055 2015 secondary 11.7240 0.0000 4.7346 736669.75
## 1 2 1202 2015 primary 27.0008 0.0000 0.0000 511589.79
## 2 3 1120 2015 secondary 11.7240 0.0000 4.7346 75562.60
## 3 4 1430 2015 secondary 0.0000 16.8225 0.0000 1498505.92
## 4 5 1041 2015 primary 11.7240 0.0000 26.2766 35298.01
## 5 6 1132 2015 secondary 11.7240 0.0000 4.7346 224773.89
## Shape_leng X_COORD Y_COORD SOURCE
## 0 7741.335 268429.9 5425167 WorldView-01 Satellite imagery from Digital Globe
## 1 5184.650 295750.8 5413701 WorldView-01 Satellite imagery from Digital Globe
## 2 1255.828 268868.3 5421731 WorldView-01 Satellite imagery from Digital Globe
## 3 15290.128 303278.4 5385710 WorldView-01 Satellite imagery from Digital Globe
## 4 2259.001 273823.9 5427266 WorldView-01 Satellite imagery from Digital Globe
## 5 4090.754 276459.4 5419663 WorldView-01 Satellite imagery from Digital Globe
## CLASSIFICA OWNERSHIP SrcOth_nad SrcOth_dat
## 0 main body of glacier Glacier National Park 0.0000 9999/01/01
## 1 main body of glacier Glacier National Park 18.6183 2014/10/19
## 2 main body of glacier Glacier National Park 0.0000 9999/01/01
## 3 main body of glacier Glacier National Park 0.0000 9999/01/01
## 4 main body of glacier Glacier National Park 0.0000 9999/01/01
## 5 main body of glacier Glacier National Park 0.0000 9999/01/01
## COMMENT Src2_0912
## 0 used both 20150822 and 20150912 images for full coverage primary
## 1 20150822 only 2015 image, but high nadir, also used 20141019 image n/a
## 2 20150822 used where shading was heavy in 20150912 image primary
## 3 used color 20150822 as secondary, some offset from 20150925 image n/a
## 4 used image from 20150822 since later image had high off nadir angle secondary
## 5 no comment primary
## Src3_0925 GLACNAME SOURCE_SCA
## 0 n/a Agassiz Glacier 1:12000
## 1 n/a Ahern Glacier 1:12000
## 2 n/a Baby Glacier 1:12000
## 3 primary Blackfoot Glacier 1:12000
## 4 n/a Boulder Glacier 1:12000
## 5 n/a Carter Glacier 1:12000
This table allows for data describing a spatial object to be stored. You can work with this table just like any other table, but you will be able to display objects in it spatially.
#polygons stores the coordinates for drawing the polygons
g2015@polygons[[1]]
You’ll notice there can be multiple polygons within a polygon. This is because there may be multiple, unconnected shapes that are considered to have the same characteristics and share the same data in the table. You’ll also see all the items needed for drawing a polygon including the coordinates, any holes to put in the polygon, and coordinates for the center. For a vector object, you can find the projection info in an object called proj4string
g1966@proj4string
## CRS arguments:
## +proj=utm +zone=12 +datum=NAD83 +units=m +no_defs
Now let’s make a map of the glaciers and view them by name.The spplot function allows you to map vector data and show different colors for a data value. Let’s look at the 1966 glaciers by glacier name.
spplot(g1966, "GLACNAME")
Before you get started, there’s one other issue. The glacier names in 2015 don’t match the rest of the shapefiles. Compare 2015 to 1966 for reference:
#check glacier names
g1966@data$GLACNAME
## [1] Agassiz Glacier Ahern Glacier Baby Glacier
## [4] Blackfoot Glacier Boulder Glacier Carter Glacier
## [7] Chaney Glacier Dixon Glacier Gem Glacier
## [10] Grant Glacier Grinnell Glacier Harris Glacier
## [13] Harrison Glacier Herbst Glacier Hudson Glacier
## [16] Ipasha Glacier Jackson Glacier Kintla Glacier
## [19] Logan Glacier Lupfer Glacier Miche Wabun Glacier
## [22] N. Swiftcurrent Glacier Old Sun Glacier Piegan Glacier
## [25] Pumpelly Glacier Rainbow Glacier Red Eagle Glacier
## [28] Salamander Glacier Sexton Glacier Shepard Glacier
## [31] Siyeh Glacier Sperry Glacier Stanton Glacier
## [34] Swiftcurrent Glacier Thunderbird Glacier Two Ocean Glacier
## [37] Vulture Glacier Weasel Collar Glacier Whitecrow Glacier
## 39 Levels: Agassiz Glacier Ahern Glacier Baby Glacier ... Whitecrow Glacier
g2015@data$GLACNAME
## [1] Agassiz Glacier Ahern Glacier
## [3] Baby Glacier Blackfoot Glacier
## [5] Boulder Glacier Carter Glacier
## [7] Chaney Glacier Dixon Glacier
## [9] Gem Glacier Grant Glacier
## [11] Grinnell Glacier Harris Glacier
## [13] Harrison Glacier Herbst Glacier
## [15] Hudson Glacier Ipasha Glacier
## [17] Jackson Glacier Kintla Glacier
## [19] Logan Glacier Lupfer Glacier
## [21] Miche Wabun North Swiftcurrent Glacier
## [23] Old Sun Glacier Piegan Glacier
## [25] Pumpelly Glacier Rainbow Glacier
## [27] Red Eagle Glacier Salamander Glacier
## [29] Sexton Glacier Shepard Glacier
## [31] Siyeh Glacier Sperry Glacier
## [33] Stanton Glacier Swiftcurrent Glacier
## [35] Thunderbird Glacier Two Ocean Glacier
## [37] Vulture Glacier Weasel Collar Glacier
## [39] Whitecrow Glacier
## 39 Levels: Agassiz Glacier Ahern Glacier Baby Glacier ... Whitecrow Glacier
You’ll notice Miche Wabun is missing Glacier in its name and North Swiftcurrent Glacier should be abbreviated as N.. Fix that so you can use the glacier names as an indexing variable across years.
#fix glacier name so that it is consistent with the entire time period
g2015@data$GLACNAME <- ifelse(g2015@data$GLACNAME == "North Swiftcurrent Glacier",
"N. Swiftcurrent Glacier",
ifelse( g2015@data$GLACNAME == "Miche Wabun",
"Miche Wabun Glacier",
as.character(g2015@data$GLACNAME)))
It’s hard to get a reference for the glaciers just from the polygon file. Let’s read in an image collected by the Landsat satellite on September 5, 2014. This will provide some reference for the location. The landsat sensors measure the light that is reflected up to the top of the atmosphere. If we overlay the red, green, and blue visible bands of light, it will look just like a photograph of the earth surface. Here you will work with the raster package to read in these .tif files. These files are different from regular photograph tif files since they are actually geotiff files that contain spatial reference information.
#read in rgb imagery from landsat
redL <- raster("data\\glacier_09_05_14\\l08_red.tif")
greenL <- raster("data\\glacier_09_05_14\\l08_green.tif")
blueL <- raster("data\\glacier_09_05_14\\l08_blue.tif")
You can check the coordinate system information in the raster files by using the following code:
#check coordinate system
redL@crs
## CRS arguments:
## +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs
You’ll notice the coordinate system is the same as the vector data above. This means that you can plot them together without having to reproject. It will be easier to plot these rasters together if we turn them into a brick. A brick is a series of raster files with the same extent and resolution. This object allows us to readily overlay the objects.
#make a brick that stacks all layers
rgbL <- brick(redL, greenL, blueL)
Now that the raster files are grouped together, make a map with both the polygons. You’ll use plotRGB to make a color image of the raster. Using the plot function adds the polygons to the map. Notice the units of the x and y axis here. When you plot polygons, you can specify a single color for a polygon or if you create a vector of colors that is the same length of the total number of polygon objects (39 here).
#plot with color
#show axes for reference
#add contrast to the imagery to see it better
par(mai=c(1,1,1,1))
plotRGB(rgbL, stretch="lin", axes=TRUE)
#add polygons to plot
plot(g1966, col="tan3", border=NA, add=TRUE)
plot(g1998, col="royalblue3", add=TRUE, border=NA)
plot(g2005, col="darkgoldenrod4", add=TRUE, border=NA)
plot(g2015, col="tomato3", add=TRUE, border=NA)
It is a little hard to see the details at the scale of the entire park. You can use the ext argument to zoom in on a specific area. ext relies on the following order: xmin, xmax, ymin, ymax.Subset the plot below to zoom in on closer on a few glaciers:
plotRGB(rgbL, ext=c(289995,310000,5371253,5400000), stretch="lin")
plot(g1966, col="palegreen2", border=NA, add=TRUE)
plot(g1998, col="royalblue3", add=TRUE, border=NA)
plot(g2005, col="darkgoldenrod4", add=TRUE, border=NA)
plot(g2015, col="tomato3", add=TRUE, border=NA)
The NDVI folder contains NDVI data collected from the MODIS sensor from 2003-2016. The NDVI data is the maximum NDVI that is observed for each year and is a derived data product generated by the USGS eModis Phenology program. The maximum NDVI will characterize the peak level of vegetation during the growing season. Let’s read in all of these files. We’ll read them in as a list and each item in our list will be one year of raster data.
#set up years to read in
ndviYear <- seq(2003,2016)
#read all files into a list
NDVIraster <- list()
for(i in 1:length(ndviYear)){
NDVIraster[[i]] <- raster(paste0("data\\NDVI\\NDVI_",ndviYear[i],".tif"))
}
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Normal Sphere (r=6370997)
## ellipsoid in CRS definition
Use the str function to take a look at what is in a single raster function. Look at the first year (2003):
str(NDVIraster[[1]])
Take a closer look at the coordinate system:
#get projection
NDVIraster[[1]]@crs
## CRS arguments:
## +proj=laea +lat_0=48 +lon_0=-113 +x_0=0 +y_0=0 +ellps=sphere +units=m
## +no_defs
This is a Lambert Azimuthal Equal Area projection that has an orgin (centers the zero point of the map at that latitude and longitude) of 48 degrees and -113 degrees. x_0 and y_0 indicate if there is a false easting or northing set (moves zero point in coordinate system, typically to keep all numbers positive). +a and +b describe the radius of the ellipsoid axes used for the datum. Units indicate the units the project is in. Many projections are in meters, but some US projections will be in feet. Note that any additional latitude arguments (lat_1, lat_2) would be associated with a standard parallel (lines of high accuracy in some projections). There are none here. This is going to be a good projection to work in since it is an equal area projection. Let’s take a closer look at the NDVI data by making a plot of it:
plot(NDVIraster[[1]])
## Vector data analysis: glacier retreat The first step to analyzing glacial retreat is to make sure the glaciers are in the right projection. The NDVI raster data is in a good projected coordinate system (PCS), and we can project the glaciers into the same PCS. The projection tool will be in the sp package for vector data (note rasters have their own projection tool in the raster package). You can just directly reference the NDVI raster projection information.
#reproject the glaciers
#use the NDVI projection
#spTransform(file to project, new coordinate system)
g1966p <- spTransform(g1966,NDVIraster[[1]]@crs)
g1998p <- spTransform(g1998,NDVIraster[[1]]@crs)
g2005p <- spTransform(g2005,NDVIraster[[1]]@crs)
g2015p <- spTransform(g2015,NDVIraster[[1]]@crs)
Note: that sometimes resizing plots with multiple sources of spatial data can have issues resizing all data layers. If it looks like the extent doesn’t match after resizing, rerun your code again after resizing.
Next let’s analyze the area of the glaciers throughout each year. There is an area column in the glacier data tables, but the metadata does not indicate how the area was measured. We will calculate the area directly for each polygon so that we know that we are consistently measuring the area of each polygon. The area function will return a vector of areas in our coordinate system units for each of the polygons.
#calculate area for all polygons
#add directly into data table for each shapefile
g1966p@data$a1966m.sq <- area(g1966p)
g1998p@data$a1998m.sq <- area(g1998p)
g2005p@data$a2005m.sq <- area(g2005p)
g2015p@data$a2015m.sq <- area(g2015p)
It will be helpful to join all of this data together into a table not associated with the shapefile so that you can analyze it. Here you will use the join from plyr.
gAllp1 <- join(g1966p@data,g1998p@data, by="GLACNAME", type="full")
gAllp2 <- join(gAllp1,g2005p@data, by="GLACNAME", type="full")
gAll <- join(gAllp2,g2015p@data, by="GLACNAME", type="full")
Let’s make a plot of the area for each glacier using this table:
plot(c(1966,1998,2005,2015),
c(gAll$a1966m.sq[1],gAll$a1998m.sq[1], gAll$a2005m.sq[1],gAll$a2015m.sq[1]),
type="b",
pch=19, col=rgb(0.5,0.5,0.5,0.5), xlim= c(1965,2016),
ylim=c(0,2000000),
ylab="Area of glacier (meters squared)",
xlab="Year")
for(i in 2:39){
points(c(1966,1998,2005,2015),
c(gAll$a1966m.sq[i],gAll$a1998m.sq[i], gAll$a2005m.sq[i],gAll$a2015m.sq[i]),
type="b",
pch=19, col=rgb(0.5,0.5,0.5,0.5))
}
It might help to visualize the glacial loss more directly. Let’s make a polygon that shows the difference in glaciers between 2015 and 1966. Here, you can use the many vector operations available in rgeos. Since we want to view the area where 1966 and 2015 don’t overlap, the gDifference function will remove any areas that overlap.
diffPoly <- gDifference(g1966p, g2015p, checkValidity = 2L)
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -79835.564606180007 107115.29284187
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -77011.733007000003 99328.148267840006
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -59187.508340740002 93811.335311250004
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -62514.363778120001 96430.893114039995
## Warning in RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid"): Ring Self-
## intersection at or near point -75367.583780639994 92136.2396056
## g2015p is invalid
## Attempting to make g2015p valid by zero-width buffering
plot(diffPoly)
What does the checkValidity argument in the gDifference function do? Why do you think we used that?
Let’s look at the difference with more context now. Let’s plot the glacier retreat with the map of NDVI showing all areas of the glacial loss in black:
#plot with NDVI
plot(NDVIraster[[13]], axes=FALSE, box=FALSE)
plot(diffPoly,col="black", border=NA,add=TRUE)
Let’s take a look at patterns in maximum NDVI from 2003-2016 to see if there is an increase in vegetation throughout the time period. Let’s start by examining how NDVI might be changing within the area that glaciers are retreating. We can get all of the values that are within a polygon using the extract function. Let’s extract each year’s NDVI in the area that experienced retreat from 1966 to 2015 and take the average. We’ll have to use a for loop to work with all values in the dataset.
#extract NDVI values
NDVIdiff <- list()
meanDiff <- numeric(0)
#loop through all NDVI years
for(i in 1:length(ndviYear)){
#get raster values in the difference polygon
NDVIdiff[[i]] <- extract(NDVIraster[[i]],diffPoly)[[1]]
#calculate the mean of the NDVI values
meanDiff[i] <- mean(NDVIdiff[[i]], na.rm=TRUE)
}
Let’s take a look at the data:
plot(ndviYear, meanDiff, type="b",
xlab= "Year",
ylab="Average NDVI (unitless)",
pch=19)
Although 2015 had very high NDVI, it doesn’t appear that there is much of a trend.This is not surprising because it may take many years for vegetation to move in and fill in a newly bare glacial area. There may also be prolonged snow cover in many of these areas even if there is not a glacier there. Let’s see if vegetation is increasing in areas surrounding the glaciers. First it will be helpful to know how much NDVI is changing in each pixel. We can fit a slope using the lm function to each cell in the pixel. We can use the calc function which applies a function to every pixel in the raster. Since you will want to apply the function to every cell in each year and for every yearly raster, you will want to be sure that R knows the raster is a stack (similar to brick) of rasters.
#designate that NDVIraster list is a stack
NDVIstack <- stack(NDVIraster)
#set up lm function to apply to every cell
#where x is the value of a cell
#need to first skip NA values (like lakes)
#if NA is missing in first raster, it is missing in all
#so we can tell R to assign an NA rather than fitting the function
timeT <- ndviYear
fun <- function(x) {
if(is.na(x[1])){
NA}else{
#fit a regression and extract a slope
lm(x ~ timeT)$coefficients[2] }}
#apply the slope function to the rasters
NDVIfit <- calc(NDVIstack,fun)
#plot the change in NDVI
plot(NDVIfit, axes=FALSE)
Let’s examine how this change might be occurring around the glaciers. We can use the gbuffer function to look at the area that is 500 meters out from the largest glacier extent to make sure we are outside of the glacier. Let’s examine if there are any individual differences in each glacier.
#buffer glaciers
glacier500m <- gBuffer(g1966p,#data to buffer
byid=TRUE,#keeps original shape id
width=500)#width in coordinate system units
We can use zonal statistics to get a statistic applied across an index, but first we need to convert our vector data to raster. The zonal function only works with two rasters.
#convert to a raster
buffRaster <- rasterize(glacier500m,#vector to convert to raster
NDVIraster[[1]], #raster to match cells and extent
field=glacier500m@data$GLACNAME, #field to convert to raster data
background=0)#background value for missing data
plot(buffRaster)
We also need to remove the actual glacier from our statistics. Let’s rasterize the largest glacial extent. Since the IDs will be the same, we can subtract the two rasters to end up with an id for only in areas with the buffer.
#rasterize gralciers
glacRaster <- rasterize(g1966p, NDVIraster[[1]], field=g1966p@data$GLACNAME, background=0)
#subtract buffer from original glacier
glacZones <- buffRaster - glacRaster
plot(glacZones)
Now you can get the statistics for the rate of vegetation change in the area around the rasterized buffer using zonal statistics:
meanChange <- zonal(NDVIfit, #NDVI function to summarize
glacZones,#raster with zones
"mean")#function to apply
head(meanChange)
## zone mean
## [1,] 0 0.0011640432
## [2,] 1 0.0008101119
## [3,] 2 0.0001584056
## [4,] 3 0.0018559251
## [5,] 4 -0.0014669432
## [6,] 5 0.0027154395